Algunos de los materiales usados en estas notas fueron facilitados por Udacity y los doctores David Luebke y John Owens, dichos materiales son parte del curso "Intro to Parallel Programming", si al finalizar de leer las notas el lector desea continuar su aprendizaje de CUDA le sugerimos ir a Udacity, además de acudir a los materiales adicionales que se encontrarán al final de cada uno de los capítulos.
Que la fuerza los acompañe.
Nuestra primer misión en CUDA es explicar ciertas sutilezas, para esto iniciaremos con un pequeño ejemplo.
Supongamos que queremos elevar al cuadrado los números del 1 al 100. Un programa en python podríamos hacerlo de la siguiente manera.
Por último lo imprimos para mostrarlo al lector.
In [1]:
enteros = range(1,101)
cuadrados = []
for i in enteros:
cuadrados.append(i*i)
print(cuadrados)
Como podemos ver, también hemos calculado el tiempo que se toma en hacer todo el proceso. Analicemos un poco nuestro algoritmo; el programa hace lo siguiente
Supongamos a manera de esquema que cada multiplicación en nuestro procesador toma 2 ns en realizarse, como estamos haciendo 100 multiplicaciones, entonces nos tomaría 200 ns obtener todos y cada uno de los cuadrados que queremos.
¿Cómo podríamos hacer esto en paralelo?
Bueno, tal vez utilizando 100 trabajadores (a los cuales de aquí en adelante llamaremos threads o hilos) y haciendo que cada uno se encargue de leer un número del arreglo de entrada, multiplicarlo por sí mismo y escribirlo dentro de otro arreglo. Si vemos lo que hace un sólo thread no podríamos notar cuál es la diferencia entre el algoritmo en Python (algoritmo escrito en serie) y el algoritmo en parelelo. La principal diferencia radica en que en el algoritmo serial un solo thread se encarga de multiplicar todos los números, avanzando de uno en uno, mientras que en el algoritmo en paralelo le damos un número a cada uno de los 100 trabajadores y los hacemos multiplicarlos por sí mismos. Si cada uno de estos trabajadores tardara 10 ns en hacer una multiplicación ¿Cuánto se tardarían en multiplicar los 100 números? Exacto, 10 ns.
Ahora veamos cómo quedaría éste código escrito en CUDA, recordando que al final de cada declaracion se debe poner un punto y coma.
#include <stdio.h>
__global__ void elevar_al_cuadrado(float * d_salida, float * d_entrada){
int idx = threadIdx.x;
float f = d_entrada[idx];
d_salida[idx] = f*f;
}
int main(int argc, char ** argv) {
const int TAMANIO_ARREGLO = 100;
const int BYTES_ARREGLO = TAMANIO_ARREGLO * sizeof(float);
// Generamos el arreglo de entrada en el anfitrion
float h_entrada[TAMANIO_ARREGLO];
for (int i = 0; i < TAMANIO_ARREGLO; i++) {
h_entrada[i] = float(i);
}
float h_salida[TAMANIO_ARREGLO];
// Declaramos apuntadores de memoria en GPU
float * d_entrada;
float * d_salida;
// Reservamos memoria del GPU
cudaMalloc((void**) &d_entrada, BYTES_ARREGLO);
cudaMalloc((void**) &d_salida, BYTES_ARREGLO);
// Copiamos informacion al GPU
cudaMemcpy(d_entrada, h_entrada, BYTES_ARREGLO, cudaMemcpyHostToDevice);
// Lanzamos el kernel
elevar_al_cuadrado<<<1, TAMANIO_ARREGLO>>>(d_salida, d_entrada);
// Copiamos el arreglo resultante al GPU
cudaMemcpy(h_salida, d_salida, BYTES_ARREGLO, cudaMemcpyDeviceToHost);
// Imprimimos el arreglo resultante
for (int i =0; i < TAMANIO_ARREGLO; i++) {
printf("%f", h_salida[i]);
printf(((i % 4) != 3) ? "\t" : "\n");
}
cudaFree(d_entrada);
cudaFree(d_salida);
return 0;
}
Para compilar este programa lo guardamos en un archivo con extensión .cu, por ejemplo cuadrado.cu, abrimos una terminal en la carpeta en que se encuentra el archivo y tecleamos el siguiente comando (el signo ! se utiliza en el notebook para hacer uso de comandos de shell, es decir, mandar comandos a la terminal desde aquí)
Nb. Guardaremos los programas en la carpeta llamada del mismo modo
In [1]:
!nvcc -o ./Programas/cuadrado ./Programas/cuadrado.cu
Esto nos creará un objeto ejecutable de nombre cuadrado, en este campo puede ir el nombre que se deseé, no es necesario que lleve el mismo nombre que el archivo que contiene al programa.
Para ejecutarlo, abrimos una terminal en la carpeta donde se creó el objeto y tecleamos
In [2]:
!./Programas/cuadrado
No tan simple como en Python ¿cierto?
La parte
#include <stdio.h>
simplemente agrega a nuestro programa la librería básica de C, que contiene funciones como printf(), scanf(), etc. Esto no es en el sentido extricto algo relacionado con CUDA, sin embargo, al estar basado en C CUDA sigue reglas similares.
Lo siguiente en nuestro código es
__global__ void elevar_al_cuadrado(float * d_salida, float * d_entrada){
int idx = threadIdx.x;
float f = d_entrada[idx];
d_salida[idx] = f*f;
}
Esto es lo que conocemos como kernel; en el kernel se ponen los comandos que el programador desea que ejecute cada uno de los threads (nuestros Oompa Loompas del GPU). A primera vista, y como ya habíamos adelantado, las órdenes parecen ser de un código serial, en ninguna parte podemos ver exactamente donde yace el paralelismo, pero hay que tener calma, llegaremos a ese punto en este mismo notebook.
Vayamos paso a paso, esta parte del código simplemente define una función para elevar al cuadrado números en un arreglo. El lector familiarizado con C habrá notado una peculiaridad (y el lector no familiarizado está a punto de enterarse), en C no existe tal cosa como __global__ (con dos guiones bajos a cada lado); nos encontramos pues ante la primer característica de CUDA.
Al ser un lenguaje de programación diseñado para controlar tanto al CPU como al GPU CUDA debe saber de alguna forma si la función que estamos escribiendo se va a ejecutar en uno u otro. Es por eso que cuenta con ciertos calificadores de tipo de función. Los tres principales son:
__global__ : especifica que la función es un kernel. Dicha función es
Cuenta también con ciertas cualidades; debe ser una función de tipo void, es decir, que no regresa nada. Y además las invocaciones a __global__ son asíncronas, en otras palabras, la función devuelve el control al programa principal antes de haber terminado de hacer el trabajo.
__device__ : este calificador declara una función que es
__host__ : declara una función
Lo siguiente en la lista es void, esto lo que quiere decir es que la función no va a regresar ningún tipo de dato, o lo que es lo mismo, sólo va a trabajar. En ocasiones es fácil confundirse con este concepto, si no va a regresar nada ¿por qué hacerlo?. Pongamos dos situaciones distintas, en la primera quieres saber cuál es el resultado de hacer cierta operación e imprimirlo, por ejemplo saber cuánto es 4!. Pero bien puede suceder que después necesites ese valor. En el caso en que sólo deseas saber cuánto vale, el tipo de función que se necesita es void, ya que sólo queremos que calcule 4! y nos diga cuánto es. En el segundo caso usaríamos una función de tipo int (integer) que regrese un número entero. La diferencia sería que la segunda función puede ser utilizada como cualquier número entero, la podemos multiplicar, sumar, dividir... A la primera no.
Después tenemos el nombre de la función, que en nuestro caso es elevar_al_cuadrado y recibe dos argumentos, declarados como
float * d_entrada
float * d_salida
Por convención, las variables que serán usadas por el dispositivo (device) se marcan con el prefijo d_
mientras que las utilizadas por el anfitrión (host) llevan el prefijo h_. Así pues, vemos que los dos argumentos que recibe nuestra función son bichos alojados en el dispositivo, dos bichos que reconocemos como números de punto flotante, pero; hay un asterisco que nos estorba, ¿será esto causante de una nueva sección en las notas?
Un apuntador es una variable cuyo valor es la dirección en memoria de otra variable. Supongamos que tenemos lo siguiente en C
float pi = 3.1416;
float * apuntador_pi;
apuntador_pi = π
En el primer renglón declaramos un flotante llamado pi con cierto valor, en el siguiente declaramos un apuntador que será dirigido hacia un número flotante, y por último dirigimos el apuntador hacia el valor deseado, usando el caracter &.
Si quisiésemos imprimir el valor de apuntador_pi la computadora escupiría algo como bffd8b3c mientras que al imprimir *apuntador_pi obtenemos 3.1416.
La pregunta que nace naturalmente es, ¿por qué complicarse la vida con apuntadores? Bueno, los apuntadores son necesarios cuando se necesita reservar memoria dinámicamente, o cuando se manejan cantidades grandes de información. Los apuntadores fueron creados en una época en que se necesitaba un uso de memoria rápido y eficiente dado que las máquinas eran demasiado lentas, eso quiere decir que en la actualidad podemos usar apuntadores para hacer correr más rápido nuestros programas o usar menos memoria, esto último, como veremos más adelante es una característica deseable cuando utilizamos GPU.
Entonces, hasta ahora hemos visto que en el primer renglón se declaró una función, ejecutable en el dispositivo pero invocable solamente desde el anfitrión, que no regresa ningún tipo de dato, llamada elevar_al_cuadrado, que recibe como argumentos dos apuntadores a estructuras alojadas en el dispositivo (demasiado para un renglón si me preguntan).
Lo siguiente en la lista es
int idx = threadIdx.x;
Declaramos una variable llamada idx y le asignamos el valor de threadIdx.x, como recordará el lector, a nuestros trajadores le habíamos cambiado el nombre por threads, threadIdx.x simplemente nos dice cuál es la etiqueta que porta el thread, o el número de trabajador si se quiere ver de otra forma; dicho número es exclusivo para cada uno de los threads, y como podemos ver, con el uso del sufijo .x los threads tienen etiquetas en cada dimensión. Es decir que podríamos acomodar threads en una cuadrícula tridimensional en donde cada uno estaría al tanto de su posición en x, y, z.
Finalmente tenemos
float f = d_entrada[idx];
d_salida[idx] = f*f;
En la primer parte declaramos una variable f que guardará el valor de la entrada ubicada en la posición idx del arreglo cuyo apuntador se metió como argumento a la función. En la segunda elevamos el número al cuadrado y lo guardamos en la entrada idx del segundo arreglo cuyo apuntador se ingresa en la funcion como argumento.
Como prometimos esta función (o kernel), será ejecutado en el GPU. Le pedimos (de una forma conveniente) 100 trabajadores al dispositivo y hacemos que cada uno de ellos ejecute el kernel. Ya antes mencionamos que cada uno de ellos contará con un índice (número de trabajador) único, y por lo tanto leerán una entrada específica del arreglo de entrada y escribirán en una entrada específica del arreglo de salida, por lo tanto no habrá que temer que dos threads intenten leer y escribir en los mismos sitios. Más adelante veremos que éste tipo de proceso es conocido como mapeo en el contexto de patrones de comunicación.
Hasta ahora hemos visto cómo se especifica cuál es el trabajo que hará cada uno de los threads, sin embargo, en ningún momento los pusimos a trabajar ¿cómo es que se logra esto? ¿nuestro kernel sirve para el número de threads que deseemos? Es hora de explicar con detenimiento el proceso de ejecución del kernel.
Primero declaramos la función main (función principal) con dos argumentos, argc y argv
int main(int argc, char ** argv)
Como vimos anteriormente, para lanzar un programa se utiliza una instrucción tipo
./nombre_del_programaSin embargo los programas también pueden ser lanzados con argumentos, supongamos que queremos un programa que promedie tantos números como le pasemos en la linea de comandos, una opción (que puede resultar tediosa) es usar scanf() y printf() cuantas veces sean necesarias, la otra es pasarlos como argumentos en la linea de comandos, por ejemplo
./promediar 2 3 4 5podría funcionar, la forma en que nuestro programa sabe cuántos números va a promediar es por medio de argc (argument count), esta variable dice cuántos argumentos se le dieron al programa. Para saber qué números va a promediar se utiliza argv (argument vector), en nuestro caso por ejemplo tendríamos
argv[0] = 2
argv[1] = 3
argv[2] = 4
argv[3] = 5
Además vemos que en los argumentos aparece char ** argv, no es porque sea un súper apuntador; con el doble asterisco indicamos que le estamos pasando un apuntador a un arreglo de apuntadores.
Si el programa no requiere argumentos es equivalente a escribir main() sin embargo es una buena costumbre agregar estos campos.
La siguiente parte del código es simple de entender, más no intuitiva.
const int TAMANIO_ARREGLO = 100;
const int BYTES_ARREGLO = TAMANIO_ARREGLO * sizeof(float);
El primer renglón simplemente define el tamaño del arreglo que usaremos. En nuestro caso, dado que queremos obtener los cuadrados de los 100 primeros números naturales (consideramos que el 0 es un número natural) declaramos una constante que guarde el tamaño del arreglo. En el segundo renglón obtenemos el tamaño en bytes del arreglo y lo guardamos en una nueva variable. Es importante observar aquí que no supusimos ningún tamaño específico para los números flotantes, sino que lo obtuvimos utilizando la función sizeof(), al momento de programar es importante escribir de tal manera que nuestros programas puedan ser corridos en otras computadoras.
Quizá se estarán preguntando ¿para qué guardar el tamaño del arreglo en bytes? Puede parecer un poco inútil, sin embargo detrás de éste renglón se esconde una sutileza importante de CUDA. No, no es el hecho de que escribimos "ni" en lugar de "ñ"
Después declaramos un arreglo de números flotantes, lo llenamos con los números del 0 al 99, y declaramos otro arreglo de flotantes para almacenar en él los cuadrados.
float h_entrada[TAMANIO_ARREGLO];
for (int i = 0; i < TAMANIO_ARREGLO; i++) {
h_entrada[i] = float(i);
}
float h_salida[TAMANIO_ARREGLO];
Podría el lector decir ¿quién utilizará esos arreglos?
float * d_entrada;
float * d_salida;
cudaMalloc((void**) &d_entrada, BYTES_ARREGLO);
cudaMalloc((void**) &d_salida, BYTES_ARREGLO);
Ahora estamos nos enfrentamos a otra de las sutilezas de CUDA, los primeros dos renglones están claros, definimos dos apuntadores de tipo flotante para los arreglos en el dispositivo. Pero ¿qué son los dos renglones siguientes? Bueno, en CUDA hay que pedir la memoria que vamos a ocupar; si bien al principio puede resultar una tarea tediosa nos ayuda a lograr el manejo eficiente de memoria que deseamos, logrando con esto tener programas más rápidos y que no consuman tantos recursos de la máquina. Para pedir la memoria que necesitamos utilizamos la funcion cudaMalloc (viene de memory allocation), que recibe como argumentos para quién reservamos la memoria y cuánta memoria necesita.
Ahora aparece un doble apuntador (void**). Todas las funciones de CUDA regresan un código de error (o uno de logro si es que no hubo errores). Todos los demás parámetros son pasados por referencia. Sin embargo, en C no se puede tener referencias, ésta es la razón por la cuál se debe pasar una dirección a la variable que queremos que guarde la información deseada. Dado que estamos regresando un apuntador, necesitamos pasarle un doble apuntador. Esto por el hecho de que cudaMalloc define un apuntador al principio, si pasáramos un apuntador simple al regresar valores no estaríamos grabando en un apuntador sino en una variable, y no es lo que queremos. CUDA necesita modificar el apuntador dado, es por esto que necesitamos pasarle un apuntador que apunte al apuntador que queremos modificar. Al final del capítulo se anexan materiales para el lector que deseé adentrarse más en este tema.
Uno debe tener cuidado al tratar de manejar las variables que hayamos declarado, debido al hecho de que el CPU no puede manejar datos alojados en el GPU y viceversa, esta es la razón por la cual existe la convención de poner h_ o d_ dependiendo de dónde planeemos utilizar la variable. Sin embargo, obviamente hay forma de compartir dados entre CPU y GPU, esto mediante la función cudaMemcpy() (de memory copy), como argumentos recibe
En nuestro programa copiamos h_entrada hacia d_entrada, le pasamos el tamaño en bytes del arreglo y por último le especificamos que la copia se hará desde algo que se encuentra en el CPU hacia algo en el GPU.
cudaMemcpy(d_entrada, h_entrada, BYTES_ARREGLO, cudaMemcpyHostToDevice);
Los primeros tres campos son simples, sin embargo en el cuarto hay que hacer una aclaración; se puede copiar de donde sea hacia donde sea; esto es:
Desde este punto de vista la programación en CUDA es heterogénea, separamos de manera tajante lo que se encuentra, manipula y ejecuta en cada una de las partes.
Y así, llegamos a la parte que todos han estado esperando
elevar_al_cuadrado<<<1, TAMANIO_ARREGLO>>>(d_salida, d_entrada);
¿Qué es esta cosa misteriosa? Bueno, el hecho de que reciba argumentos nos hace pensar que es una función (y lo es), el nombre elevar_al_cuadrado al principio nos indica que nuestra corazonada es correcta, es la función que definimos al principio; nuestro kernel, esa lista misteriosa de trabajos para los threads en donde yacía el paralelismo que no nos era posible observar. Bueno, ahora podemos verlo. Además de las características que mencionamos tiene una resaltable, aparece el símbolo < tres veces después del nombre de la función ¿quiere decir esto que elevar_al_cuadrado es muchisísimo menor que 1? La respuesta es, no. Resulta que los kernels reciben un argumento extra, de alguna manera tenemos que decirles cuántos threads queremos que ejecuten el kernel, eso lo hacemos utilizando <<<1,TAMANIO_ARREGLO>>>, de esta forma le estamos diciendo, quiero que lances una fila con un número de trabajadores igual a lo guardado en la variable TAMANIO_ARREGLO (en nuestro caso 100) para que ejecuten el kernel.
La respuestas se encuentran en...
En CUDA se tienen 3 tipos principales de estructuras en lo que se refiere a lanzar kernels.
Thread: Es simplemente un ente encargado de la ejecución del kernel, con un índice dado. Cada thread conoce su índice para acceder a elementos en un arreglo, de tal manera que la colectividad de todos los threads procesa cooperativamente toda la información.
Block (bloque): Es un grupo de threads. No hay demasiado que decir acerca de la ejecución de los threads dentro de un bloque - pueden ejecutarse serialmente o con ningún orden particular. Se pueden coordinar los threads, de cierta forma, usando una función llamada _syncthreads() (que veremos más adelante) que hace a los threads detenerse en cierto punto de la ejecución del kernel y continuar hasta que todos los otros threads del bloque hayan llegado a dicho punto.
Grid (malla): Es un grupo de bloques. No existe sincronización alguna entre bloques.
Ahora, tómense fuerte del asiento porque estamos a punto de hacer un recorrido dimensional muy interesante.
Al lanzar el kernel elevar_al_cuadrado lo lanzamos con nuestra notación muchisísimo menor que <<< >>> y como argumentos le pasamos 1, TAMANIO_ARREGLO lo que hicimos aquí fue especificar que deseábamos un bloque de $1\times100$ threads, un arreglo unidimensional. Pero, ¿que tal si quisiésemos procesar una imagen? dado que las imágenes son puramente bidimensionales sería mejor lanzar un kernel de dos dimensiones, por ejemplo si tuviésemos una imagen de $256\times256$ una de las apuestas que podríamos hacer sería lanzar un bloque así <<<256,256>>>. Bueno, la diversión no termina ahí, podemos hacer bloques de tres dimensiones también <<<256,256,256>>>, ¿y las mallas? En estos casos estámos especificando tácitamente que no hay malla, o de otra forma, que nuestra malla es de $1\times1$. Esto nos lleva a lo siguiente, podemos lanzar también arreglos de mallas. La dimensión de la malla y de los bloques son variables del tipo dim3 así que por ejemplo el kernel de nuestro ejemplo pudimos haberlo lanzado así
dim3 bloque(1,TAMANIO_ARREGLO);
elevar_al_cuadrado<<<bloque>>>(d_salida,d_entrada);
De la misma manera podemos definir mallas de dos o en algunos casos 3 dimensiones y lanzarlas así
dim3 malla(4,12);
dim3 bloque(256,256,256);
kernel<<<malla,bloque>>>()
Hay dos limitantes que hay que tener en cuenta
Por último, el thread conoce exactamente sus coordenadas dentro del bloque, las coordenadas del bloque dentro de la malla, y las dimensiones tanto del bloque como de la malla. Todas estas propiedades pueden ser accedidas con las variables
threadIdx.x.blockIdx.x, blockDim.x.gridDim.x.En donde en cada caso podemos sustituir la x por y,z. Más adelante veremos que todas estas variables son útiles al momento de asigarle indices globales a cada thread.
El lector a estas altura es capaz de decir qué hace la siguiente linea de código
cudaMemcpy(h_salida, d_salida, BYTES_ARREGLO, cudaMemcpyDeviceToHost);
Pero por el placer de llenar espacio lo recordamos, se llama a la función cudaMemcpy() para que transfiera hacia el arreglo h_salida los valores de d_salida, que tiene un tamaño en bytes BYTES_ARREGLO y especificamos que la copia se hará desde el GPU hacia el CPU utilizando cudaMemcpyDeviceToHost.
Imprimimos el arreglo como se hace usualmente en C
for (int i =0; i < TAMANIO_ARREGLO; i++) {
printf("%f", h_salida[i]);
printf(((i % 4) != 3) ? "\t" : "\n");
}
Y para terminar liberamos la memoria con
cudaFree(d_entrada);
cudaFree(d_salida);
return 0;
y regresamos un bonito 0 para indicar que todo ha ido de maravilla, si no devuelve un 0 la función main entonces el programa sabrá que hubo un error.